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The free cooling process in the inelastic hard sphere system is studied by analysing the 
data from large scale molecular dynamics simulations on a three dimensional system. The ini- 
tial energy decay, the velocity distribution function, and the velocity correlation functions are 
calculated to be compared with theoretical predictions. The energy decay rate in the homoge- 
neous cooling state is slightly but distinctively smaller than that expected from the indepen- 
dent collision assumption. The form of the one particle velocity distribution is found not to 
be stationary. These contradict to the predictions of the kinetic theory based on the Enskog- 
Boltzmann equation and suggest that the velocity correlation is already important in the early 
stage of homogeneous cooling state. The energy decay rate is analysed in terms of the velocity 
correlation. 
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1. Introduction 

The hard sphere system is the oldest system studied 
in the statistical mechanics and is still attracting much 
interest of current research. Some of the recent studies 
on this system were triggered by the recent appreciation 
of the fact that, no matter how small the inelasticity is, it 
could cause drastic changes in the system behavior and 
leads to series of transitions as the system evolves, as 
long as the system is large enough.^) 

In the inelastic system, the particles lose relative veloc- 
ity every time the particles collide each other. The first 
consequence of this is the uniform cooling of the sys- 
tem (Homogeneous Cooling State, HCS). In this stage, 
no macroscopic structure is assumed to exist and we can 
measure the kinetic energy in terms of the granular tem- 
perature T defined by 



(1) 



where d is the system dimensionality, m is the mass, and 
(• • • ) represents the statistical average. Then it has been 
shown that the initial cooling follows as 



Tn 



(2) 



(l + iAo)^' 

where t is the time, and and Tq are constants. This is 
called Haff's law. 2) 

Secondly, the velocity distribution deviates or does not 
converge to the Maxwell-Boltzmann distribution. '^^ The- 
oretical analyses and numerical simulations have shown 
that the velocity distribution in an inelastic system falls 
into a distribution different from the Maxwell-Boltzmann 
within a first several collisions per particle if the system 
is uniform."*' This form of distribution has been shown 
to be a stable and stationary solution of the Enskog- 
Boltzmann equation of the kinetic theory for an inelas- 
tic system. Molecular dynamics simulations on a two di- 
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mensional system, however, have shown that this form of 
scaled distribution does not actually remain stationary 
but gradually returns back to the Maxwell-Boltzmann 
due to the velocity correlation ignored in the Enskog- 
Boltzmann approximation.^' 

This velocity correlation in the inelastic system is de- 
veloped by the reduction of relative velocity of colliding 
particles during the global cooling. Since the inelastic col- 
lision reduces kinetic energy but not momentum, random 
fluctuations in momentum over a certain spatial region 
come to stick out when the size of momentum fluctu- 
ations becomes comparable with the momentum itself. 
This is called the noise reduction instability.^' The ve- 
locity correlation may be seen as a vortex structure in 
the velocity field of the system. 

Well developed vortex structure causes shearing in flow 
between vortices, which in turn causes viscous heating. 
The pressure gets higher in the heated region, then this 
leads to particle flow due to the resulting pressure gradi- 
ent, and eventually inhomogeneity of the system in the 
density field, '^' which is called the clustering instability. 

In this paper we present a detailed analysis of numer- 
ical simulations in a three dimensional system, and give 
evidences that velocity correlations have a deep influence 
on the behavior of inelastic hard sphere system even in 
a very early stage. 

The paper is organized as follows; In §2, we introduce 
our model system and the simulation setup. In §3, we 
present our simulation result for the time dependence of 
the kinetic energy. In §4, the velocity distribution is com- 
pared to the one obtained from the theory based on the 
Enskog-Boltzmann equation using the generalized La- 
guerre polynomial expansion to show the deviation from 
the Gaussian. In §5, we compare the velocity correlation 
with the theoretical result based on the hydrodynamic 
equation with the Langevin force. In §6, the deviation 
of the energy decay rate from Haff's law is analysed in 
detail. We summerize our results in S7. 



2. Model system 

The model system wc study is a three dimensional (3- 
d) hard sphere system which undergoes inelastic colli- 
sions with a constant normal restitution coefficient e. We 
consider the mono-dispcrsc system with the diameter g 
and the mass m. The particle rotation is ignored for sim- 
plicity. Then the collision rule is given in terms of the 
velocity of the i-th particle Vi and v\ before and after 
the collision, respectively, as 

1 



tion of random uncorrelated collision with the Maxwell- 
Boltzmann velocity distribution, 7 becomes 
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-[n- {vi-Vj)\n, 
-[n ■ {vi - Vj)]n, 



where i and j are the colliding particles, and n is the unit 
vector parallel to the relative position of the colliding 
particles at the time of contact. In the following, we use 
the unit system where the particle diameter, the mass, 
and the initial average kinetic energy per particle are 
unity. 

To simulate the system, we employ the event driven 
method, in which we keep track of colliding events. The 

actual code of the simulations is based on the fast algo- 
rithm by Isobe,^-' which enables us to simulate the system 
with the particle number as many as 10^ in 3-d on our 
PC's. Most of the simulations presented here are done 
for the restitution coefficient in the range of e > 0.8 with 
the number density n = 0.4, which corresponds with the 
volume fraction = {A/3)TT{a/2)^n = 0.209. The density 
dependence is also simulated for some observables. 

The system is contained in the cube with the peri- 
odic boundary condition. The initial state is prepared 
as the equilibrium state by running the system long 
enough with the restitution coefficient e = 1. To check 
the equilibration, we monitor the velocity distribution 
and wait until some of the Laguerre coefficients become 
zero within a statistical error; this procedure typically 
takes about 50 collisions per particle. 

The inelastic collapse'^' does not occur without any 
special treatment in the present simulations because we 
mainly focus on the initial and intermediate stage of cool- 
ing with the parameters of relatively low dissipation. 

3. Energy decay in HCS 

As the system cools, things slow down and colliding 
events become less and less frequent. Therefore, it is con- 
venient to measure the time in terms of the average num- 
ber of collisions r that each particle experiences. This is 
related with the time t by 

dr = ojdt, (3) 

where iv is the collision frequency, which depends on 
time. 

In HCS, collisions are thought to be uncorrelated and 
the average energy loss in each collision is proportional to 
the kinetic energy, or the granular temperature, therefore 
we have 

f = -2,.T, (4) 
where 7 is the energy loss rate. Within the approxima- 
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(5) 



because only the normal component of the velocity re- 
duces by the factor e. This gives Half's law for Maxwell- 
Boltzmann velocity distribution, 



THaff(T) = Toexp[-27or]. 



(6) 
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Fig. 1. The time dependence of the granular temperature T de- 
fined in eq.(l) for various restitution coefficients e with n = 0.40 
and N = dx 10^ . The average number of collisions t is scaled as 
270T . The solid straight line represents the Haif's law with 70 
(eq.(6)). 



Figure 1 shows the energy decay; the granular tem- 
perature T{t) defined in eq.(l), which is proportional to 
the kinetic energy, is plotted as a function of 270T for 
several values of e. Although the initial exponential de- 
cay seems to be well represented by Haff's law (6), there 
is a small but distinctive discrepancy in the decay rate 
7 from eq.(5) as we will see in Fig. 11 in §6, where we 
plot numerically obtained decay rate 7 with the theoret- 
ical analysis that we will discuss in detail. There arc two 
possible causes for this: (i) the deviations of velocity dis- 
tribution from the Maxwell-Boltzmann distribution, (ii) 
the correlation in collisions. We will discuss (i) and (ii) 
in §6. 

The energy decay deviates from the exponential as the 

time elapses, after then the clustering occurs. We define 
the crossover time Tc between Haff's region and the clus- 
tering region as the time when Tuasir) becomes half of 
T(r), namely 



-T(Te) = THaff(rc), 



(7) 



where Tuas{T) is given by eq.(6). In the case of n = 
0.40, we obtain = 58, 82, 134, and 309 for e = 
0.80, 0.85, 0.90, and 0.95 respectively. The crossover time 
Tc is plotted against 270 in Fig. 2. In the following sec- 
tions, we will focus on the system behavior in early stage, 

T < Tc. 

The beginning of the clustering has been observed 
around r ^ Tc in our simulation and it has a strong 
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Fig. 2. The crossover time Tc as a function of 270 = (1 
e^)/d. The circles represent our simulation results for e 
0.95, 0.90, 0.85, and 0.80 from left to right. 



effect on the energy decay. The long time behavior of the 
energy in the clustering region seems to be some power 
decay like t^". From the hydrodynamical approxima- 
tion, a = d/2 has been obtained/°^ but our result seems 
to show = 3-^4 with some sample dependence and 
parameter dependence. 

4. Velocity distribution 

The velocity distribution function /(u, t) in the inelas- 
tic system is conveniently represented in the scaled form 
using the average speed 



= \j dvf{v,T)i 



(8) 



as 



1 V 

The scaled velocity distribution function p{c, r) would be 
independent of r if the form of the velocity distribution 
were stationary as the system cools down. 

When the distribution function is not far from the 
Gaussian, it is reasonable to expand it as 

Pic r) = exp(-c2) J2 cie{r)Lnc% (a = - - 1), (10) 

e=o ^ 

using the generalized Laguerre polynomials, or the So- 
nine polynomials,^' L'^{x). The coefHcients a£{T) can 
be obtained by 



aeir) 



£lT{d/2) 



J dcp{c,T)L'^{^), 



(11) 



r(^ + d/2) 

where T{x) is the gamma function. The values of the first 
two coefficients are given by 

0,0 = 1, fli = 0, (12) 

due to the normalization and the scaling, but any devia- 
tion from the Maxwell-Boltzmann distribution is repre- 
sented by non-zero values of the higher order coefficients 

at {£ > 2). 

The time development of ae{T) has been calculated 
based on the Enskog-Boltzmann equation with a certain 



truncation scheme for the higher order coefficients. It 
has been shown that the Enskog-Boltzmann equation has 
the solution which starts from ao = 1 and = (i > 
1) and converges to non-zero constants within a several 
collisions per particle. 
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Fig. 3. The scaled velocity distribution function p(c) at a several 
values of r for the system with e = 0.80, n = 0.40, N = 6 x 
10^ . The crossover time Tc is 58 for this system. The solid line 
represents the Maxwell-Boltzmann distribution p(c) = exp(— c^). 
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Fig. 4. The time development of the generalized Laguerre poly- 
nomial expansion coefficient of the p(c) for e = 0.90, n = 
0.40, AT = 6 X 10^. The data are averaged over three realiza- 
tions. The crossover time Tc is also shown by the arrow. 



Figure 3 shows the scaled distribution of normalized 
particle speed c for the system with e = 0.8 for a several 
values of the time r in the semilogarithmic scale. The 
deviation from the Gaussian is very difficult to see even 
after the clustering sets in around Tc = 58. 

The Laguerre coefficients a2, 03, and 04 are shown as 
functions of r in Fig. 4 for e — 0.90 and n = 0.40. They 
are small but clearly have non zero values, which shows 
that the coefficients deviate from those for the Gaussian. 
The initial deviation from the Gaussian is very quick as 
has been predicted, and \a2\ > Ids] > \a4\ holds, which 
supports the validity of the truncation scheme. In Fig. 
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Fig. 5. Comparison of 02 of the simulation result (solid line) Fig. 7. Long time behavior of 02 for various values of the resti- 



witli the solution of the Enskog-Boltzmann equation with the 
second order approximation given in cq. (A-4) (dashed line) for 
e = 0.90, n = 0.40 and N = 6 X 10^ (averaged over 3 realiza- 
tions). The inset shows the short time behavior (averaged over 
6 realizations). 



tution coefficient e at n = 0.40 and N = 6 X 10''. The 
data presented in the graph are smoothed by the time aver- 
age of one sample over the interval At = 3, 5, and 7 for 
e = 0.80, 0.85, 0.90, and 0.95, respectively, from the original 
data being measured at every integer t. 



5, the simulation results are compared with the solution 
of the Enskog-Boltzmann equation for 02 (t), whose ex- 
pression is given in Appendix A. The simulation results 
agree very well with the solution till it becomes close 
to the stationary value, but after then a2(T) and other 
coefficients go back towards zero until r ~ Tc. 
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Fig. 6. The time dependence of 02 for various values of the den- 
sity n with e = 0.85 and N = 6 x 10®. The solid line represents 
the solution of the Enskog-Boltzmann equation with second or- 
der approximation given in eq.(A-4). 



The coefficients obtained from the simulations depend 
on the density although the Enskog-Boltzmann equation 
results do not. Figure 6 shows the short time behavior of 
02 (r) for various densities n. 02 (r) for the smaller den- 
sity case is close to the theoretical solution. This suggests 
that this behavior of returning to the Gaussian is caused 
by velocity correlation because multiple collisions of the 
same pairs occur more frequently in higher density. Ac- 
tually, it has been shown in the 2-d simulation that the 
Laguerre coefficients stay stationary if the velocity cor- 
relation is destroyed artificially by shuffling the velocity 
of each particle. 



The behaviors of 02 (r) for various restitution coeffi- 
cients e are shown in Fig. 7 with the time being scaled 
by the crossover time Tc- All of them show similar trend 
when the time is scaled with Tc- The behavior after the 
clustering (t > r^) has large sample dependence. This 
seems to be due to the fact that the system is not large 
enough compared with the cluster size. 

5. Velocity field and its correlation 

The random velocity field in the initial state becomes 
organized during HCS through inelastic collisions. This 
has been clearly seen in 2-d system as the vortex struc- 
ture.^^' In 3-d system, it is rather difficult to see directly 
due to a problem of visualization and also due to the 
difficulty to have a large system which contains enough 
vortices to see their structure. 

Figure Sshows a snapshot of a velocity field just be- 
fore the clustering starts. The velocity field u{r,T) is 
obtained by averaging over a velocity of all the particles 
in a coarse grained cell at r; 

-j^ cell at r 

u{r,T) = — J2 (13) 

^ i 

where Nr is the number of the particles in the cell at 
r. In this figure, the volume of the cell is (1/32)^ of the 
whole system. In order to see the 3-d flowing structure, 
the velocity field is shown only in three thin layers in 
the system. The field is rather complex, but there can be 
seen some organized vortex like structure of the hydro- 
dynamic scale. 

In order to quantify the structure in the velocity field, 
we measure the velocity correlation function defined as 

Ga/3{ri,r2; t) = (wa(''i, T)u^(r2, r)). (14) 

In the case of the isotropic uniform system, this can 
be decomposed into two parts as 

Ga0{r;T) = faf0G\\{r,T) + {da/s - far/3)G±{r,T), (15) 




Fig. 8. Cross sections of velocity field at r = 60 (tc = 58) for 
n = 0.40, e = 0.80 and AT = 1 x 10'^ with the system size L = 136. 
Only 1/4 of the xy plane and 3/32 of the z direction of the whole 
system are shown. The presented velocity field is the average 
velocity over all particles in a cell whose size is 1/32^ of the 
system. The tone and the length of the arrows represent the 
magnitude of the velocity (the darker arrow represents larger 
velocity). 
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Fig. 9. The velocity correlation Gy and Gx at t = 50 for n = 
0.40, e = 0.80 and = 1 X 10*' just before the crossover time Tc = 
58. The theoretical curves are calculated for the incompressible 
case. 




Fig. 10. Scaling form of the velocity correlation Gy (upper curves 
with monotonic decay) and Gx (lower curves with negative min- 
imums) for n = 0.40, e = 0.80 and W = 1 x 10®. 



where 

T 

r = n - r2, f = -, (16) 
r 

and the longitudinal part G\\ and the transverse part G± 
may be calculated by 

G||(r,T) = ((M(r,T)-f)(M(0,T).f)), (17) 

'^^(^'^) = rf3T<("(^'^) ■ - G||(r,T)). (18) 

These; correlations have been analyzed based on the 
hydrodynamic equation with the Langevin force. ^^^^^^ 
The actual expressions are complicated and given in Ap- 
pendix B for the 3-d case using the same notation with 
ref.ii) 

Figure 9 shows G|| and G±_ with the theoretical curves 
for the incompressible case. The agreement is quite good. 
The negative correlation in G±_ is a sign of the vortex 
structure. We also evaluated the theoretical expressions 
of G|| and G±^ for the compressible case, but the differ- 
ence between the compressible and incompressible cases 
are small for these parameters. 



The correlations G\\ and G_l for several t's are given 
in Fig. 10 using the scaling Gr'^/'^m/To v.s. r/a^Jr. All 
the plots for different r's collapse on a single curve for 
20 < T < 60, which means that the correlations scale as 

G{r,T) = ^^f(^^, (19) 

very well. We did not observe large density fluctuation 
until T - 70 (tc = 58). 

6. Correction for the energy decay rate 

As we have pointed out in §3, the energy decay rate in 
HCS slightly deviates from the one given by Haff's law. 
We shall examine this deviation in detail in this section. 

Haff's decay rate 70 is obtained on the random col- 
lision assumption with the Gaussian distribution of the 
particle velocity. We will first examine the effect of the 
non-Gaussian velocity distribution and then the velocity 
correlation. 



6.1 Effect of non- Gaussian velocity distribution 0.865. 
The time dependent decay rate defined as 

27(t) = In T(r), (20) 

and the colhsion rate uj{t) have been studied in the case 
of non- Gaussian velocity distribution using the Enskog- 
Boltzmann theory to give'^^ 



7(t) 



70 



1 + Y6«^(r) 



^o(t) 



ujo{t) 16 



(21) 



(22) 



up to the first order of 02. Here, 02 is given in Appendix 

A and luq{t) is the collision frequency in the Enskog- 
Boltzmann approximation in the equilibrium with the 
average velocity vq{t), 



'^0(Tj = 7= Vo{t), 



27r 



(23) 



where fl^ is the surface area of the d dimensional unit 
sphere and x('^) is the pair correlation at contact in the 
equilibrium with the density n. 



1 

0.99 
0.98 
0.97 
0.96 







□ □ ^ 
















E-B 
n=0.10 
n=0.20 

n=0.40 

1 


□ 


A 


■"^^^4 A 

1 



10 



15 20 



Fig. 12. Comparison of the energy decay rate 7(r) of the sim- 
ulation result (points) with theoretical estimates for various 
values of the density n with e = 0.85 and N = 0.6 x 10®. 
Density dependence is shown. The dashed hne represents the 
Enskog-Boltzmann result eq.(21) with the analytical solution of 
0,2 (eq.(A-4)). Solid lines show the phenomenological approxima- 
tion of eq.(32) with fitting parameter b = 0.30. 
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Fig. 11. Comparison of the energy decay rate 7(t) of the simula- 
tion results (points) with theoretical estimation for various values 
of the restitution coefficient e with n = 0.40 and = 0.6 X 10^. 
The dashed line represents the Enskog-Boltzmann result eq. (21) 
with the analytical solution of 02 (eq.(A-4)) for e = 0.865. Solid 
lines show the phenomenological approximation of eq.(32) with 
fitting parameter b = 0.30. 



Numerically obtained decay rates 7(t) in the early 
stage of HCS are plotted in Fig. 11 for a several values 
of restitution coefficient e. The dashed line in Fig. 11 is 
the Enskog-Boltzmann prediction (21) for e = 0.865, for 
which we expect that the deviation of 7 from 70 should 
be largest in the range of 0.80 < e < 1 because the sta- 
tionary value of \a2\ reaches the maximum at that value. 

There are two points that we can see in Fig. 11; (i) 
actual deviations of 7 from 70 are much larger than this 
largest possible deviation within the Enskog-Boltzmann 
theory, and (ii) the deviation increases monotonically 
upon increasing dissipation while the Enskog-Boltzmann 
theory expects the maximum deviation around e = 



Density dependence of 7(t) is shown in Fig. 12 for e = 
0.85. Although the stationary value of 02 does not depend 
on the density n in the theory based on the Enskog- 
Boltzmann equation, the simulation result of 7(t) docs, 
as well as the simulation results of a2- The absolute value 
of 02, however, becomes smaller than that of the Enskog- 
Boltzmann theory for larger density, while the deviation 
of the decay rate 7 from 70 gets larger, as n increases. 

These points suggest that the deviation in 7 cannot be 
explained by the non-Gaussian distribution, and should 
be mainly caused by the velocity correlation even in the 
early stage of HCS. 

6.2 Effect of velocity correlation 

Now, we will try to include the effect of velocity corre- 
lation. The velocity correlation results in the local aver- 
age flow of particles and should lead to the following two 
effects that slow down the energy decay: (i) reduction 
of collisional dissipation and (ii) viscous heating due to 
shearing of the average flow. We will examine only the 
first one since the viscous heating is the secondary effect 
and should be negligible in the early stage of HCS that 
we are looking at. 

In the presence of the local average flow u{r,T), the 
granular temperature Ti,{r,T) should be redefined as 

cell at r 

^l(^'^)^]v; E (24) 

i 

where the summation is taken over the particles in a 
cell at r. We have to use this definition instead of the 
granular temperature T given in eq.(l), which actually 
corresponds to the kinetic energy. 

Our physical assumption is that the collisional dissipa- 
tion should be determined by the granular temperature 
Tl defined in eq.(24) because the average flow cannot be 



relaxed through local collisions. Then, Half's law for the 
energy decay is extended by modifying eq.(4) as 



(25) 



where E{t) = — {d/2)T{T) is tlic average en- 

ergy per particle. 

The energy may be expressed in terms of the velocity 
correlation and the granular temperature as in the fol- 
lowing. In the presence of the average flow, the kinetic 
energy per particle E{t) can be decomposed into two 
parts as^°) 



E{t) = ^ J dr |i(mn(r,r)u2(r, 



r)} 



+ -{n{r,T)Ti^ir,T)) 



(26) 



where n(r, r) is the local particle density. 

We assume, at this point, that n and (Tl) are uniform 
because we are considering early stage of HCS, then we 
have 

Eir) « 1 y dr|(u2(r.,r)) + ^{nir)). (27) 

The first term of RHS can be expressed in terms of the 
velocity correlation in the hydrodynamic limit as 



1 1 dr^{u\r,r)) = 1 ^ lim^Tr [G(n2,T)] , 



^iTr[G+(0,r)]. 



(28) 



The correlation function has been obtained within the 
linear approximation^^' as 

rHaff(T) 



Tr[G+(0,T)] =(d-l): 



mn 



dk e2^«"(i-«i'=') - 1 



(29) 



(27r)" 
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for the incompressible case. The short wave length cutoff 
fcM is needed for the integral to give a finite result; the 
integral diverges logarithmically for 2-d and linearly for 
3-d as fcM ^ 00. We set fcM as 

fcM = 2?™^/'' X b, (30) 

where we take 6 as a fitting parameter independent of 
the density and the restitution coefficient. Then we have 

THaff(T)f2d 



Tr[G+(0,T)] =(d-l) 



mn{2'jT^±)'^ 



2-/0 r 



(31) 



ds 



2s'^/' 



n{d/2,.siiki,), 



where 7(0,, x) is the incomplete gamma function defined 
in Appendix B. For large r, this expression gives E ~ 
T"*^/^, as has been studied in the ref."'^"-' 

In the ref,^°' (Tl) is assumed to obey Half's law, but 
here we do not use this assumption; eliminating (Tl) 
from eq.(25) using eq.(27), a differential equation for 



E{t) is obtained as 

dE 



dr 



(r)=-27o(i?(r)-^Tr[G+(0,r 



(32) 



We solve eq.(32) numerically with numerically evalu- 
ated Tr [G+(0, r)] from eq.(31). 

Comparison of the energy decay rate 7(t) with the 
simulation results is shown in Figs. 11 and 12. The value 
of the fitting parameter b is chosen to be 6 = 0.30 to give 
the best fit for the energy E{t) before Tc at e = 0.80 and 
n = 0.40. Both the restitution coefficient and the density 
dependence are well described in this approximation with 
the same value of the single fitting parameter b. 
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Fig. 13. The deviation of the energy E{t) from Haff's law. The 
simulation results (points) are compared with the phenomeno- 
logical approximation of eq.(32) (solid lines) for n = 0.40 and 
AT = 0.6 X 10^. Fitting parameter is & = 0.30. The dashed line 
is HaiT's law exp(— 270T) (the crossing points of the simulation 
results and the dashed line give Tc's in our definition in eq.(7)). 



In Fig. 13, the deviation of the energy from Haff's law 
E{t) — Eiia.s{T) is plotted and the simulation results are 
compared with the results from eq.(32). The agreement 
is quite good and this approximation seems to be valid 
until little after Tc- But in the clustering region, the sim- 
ulation results have stronger energy decay than our es- 
timates. These results show that the deviation of the 
energy from Haff's law at the initial stage is mainly due 
to the velocity correlation, while the effect of the viscous 
heating and density cluster might have to be considered 
at the clustering regime. 

7. Summary and Discussions 

Unlike in the equilibrium state, the velocity correlation 
develops in the non-equilibrium state of the inelastic sys- 
tem, and it seems to play an important role in the time 
development of the system. We carefully looked into the 
early stage of the free cooling in the 3-d inelastic hard 
sphere system in the low dissipation range e > 0.8, and 
found (i) the initial energy decays exponentially, but the 
decay rate is slightly smaller than that expected in the 
case of the random collision, (ii) the velocity distribu- 
tion deviates from the Gaussian but only slightly, (iii) 
the form of the velocity distribution does not stay sta- 
tionary even in HCS contrary to the expectation by the 



Enskog-Boltzmann equation, (iv) the velocity field devel- 
ops global structure and the velocity correlations fit to 
the theoretical results based on the linear analysis of the 
hydrodynamic equation with the Langevin force, (v) the 
deviation of the energy decay rate can be described in 
terms of the velocity correlation. 

Among these results, (i), (iii), and (v) are the evidences 
that the velocity correlation plays an important role even 
in the early stage of HCS. The velocity correlation in 
this stage has also been studied numerically using the 
ring approximation of the kinetic theory, where it has 
been shown that the velocity correlation arises within a 
small number of collisions. 

The effects of the non-Gaussian velocity distribution 
and the velocity correlation on the energy decay have 
been also studied in the kinetic theory by Poschel et al}'^^ 
They studied the energy decay in the real time, and de- 
rived the expression for the deviation from Half's law 
by the non-Gaussian distribution and the velocity cor- 
relation. By comparing it with the 2-d simulation, they 
concluded that a major part of the deviation should come 
from the velocity correlation, but they did not make the- 
oretical estimate for the size of the velocity correlation. 

The effect of the velocity correlation should depend on 
the particle density; the lower the density is, the smaller 
the correlation effect becomes, because the mean free 
path gets long much faster than the mean distance be- 
tween particles upon decreasing n. However, this does 
not necessarily mean that the correlation effects disap- 
pear in HCS for small enough density; for large enough 
system, it is believed that the clustering instability takes 
places eventually for any density and dissipation. This 
should be caused by the shearing flow, and the shearing 
is a manifestation of the velocity correlation. Therefore, 
it is plausible that, even for a small density system, the 
correlation effect takes place from the early stage of HCS 
in comparison with the total period of HCS, which be- 
comes long for small n. 

In the present work, the simulations are done on 3-d 
systems and the analysis are based on the 3-d expres- 
sions, but the above findings are not particular to 3-d, 
and the physical processes do not seem to be very differ- 
ent from 2-d system, especially in the early stage where 
linear theories hold. As for the late stage cluster growth, 
some differences between 2-d and 3-d system have been 
reported from MD simulations.^''-' 

Then shouldn't there be any differences between 2-d 
and 3-d cases before the clustering? 

In the 2-d system, the velocity correlation developed 
in the system has been visualized as a vortex structure 
in the intermediate stage of HCS. The appearance of the 
structure in 2-d reminds us of the vortex structure in the 
free decaying turbulence, where the dual cascade of 
energy and enstrophy is a leading process in the inertial 
range. In the 3-d system studied here, we have not suc- 
ceeded in visualizing the clear vortex structure although 
there seems to be some global structure. If the dynam- 
ics of vortices controls the system behavior in the later 
stage of HCS, or the shearing state, there may be an 
important differences between the two cases; in 2-d, the 
vorticity a> = V x ti has only one component and the 



coarsening of the scaler field is a major process, but in 
3-d, the vortices are lines in three dimensional structure 
and their motion such as stretching and bending accom- 
panied by reconnection of vortex lines could cause more 
complex dynamics that keeps certain fine structures in 
the late stage; this is not taken in the linear theory. Such 
dynamics should have some effects in the early stage of 
clustering because the clustering is known to be triggered 
by the viscous heating in the vortex structure. 

Appendix A: Second order approximation for 
Laguerre expansion coefficients in 
higher dimensions 

In this appendix, we present the expression of a2 
within the Enskog-Boltzmann equation for the second 
order^^ for d dimensions. a2 follows the equation, 



^{r)=A + a2{T)B, 



(A-1) 



with 



A = 



(e^ - l)(2e^ - 1) 
d{d + 2) 



(A-2) 



B 



(1 + e){9 + 2Ad + 8ed - Ale + 30(1 - e)e'^) 



16d{d + 2) 



Its solution is 



a2(r) = - 
and stationary value is 



(12 



(1 - exp(Sr)), 



A 



(A-3) 
(A-4) 

(A-5) 



Appendix B: Analytical form of velocity corre- 
lation function 

In the ref.,^^) the velocity correlation functions have 
been obtained from the hydrodynamic equations with the 
linear approximation for both incompressible and com- 
pressible flows in the form of the Fourier integrals, but 
the explicit expressions for d dimensions are given only 
for the incompressible case. The expressions for d dimen- 
sions used in the text are given by 



G+ = ^ I ds'e^' 



O-A- 



{d-l^w^ 



44> 



(B-l) 



for the incompressible fluid and 



Tnaff 



270 r 



ds'e' 



{d-Vf^ 



,d X 



'^''2' 4s''' 



2' 4s' 



(B.2) 



+ 



(47rs')''/2 



for the compressible fluid, where A = {||, _L}, Sw is Kro- 
necker's delta, a = 5\\\ — 5_\_\, THaflf is the temperature 
given by Half's law (eq.(6)), 

7(a,a;)= / exp(-t)t""\ (B-3) 
Jo 

is the incomplete gamma function, 

XX = r/a , (B-4) 

and ^a's arc the correlation lengths of the elastic hard 
sphere system, given in the ref.^^-* 
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